week 2: linear model and causal inference

geocentric models

Annoucements and such

  • We will start using brms today! Install this package now, if you haven’t already.
install.packages(c("brms","tidybayes"))

Workspace setup

library(here)
library(tidyverse)
library(cowplot)
library(brms)
library(tidybayes)
library(patchwork)

Workflow

  1. State a clear question.

  2. Sketch your causal assumptions.

  3. Use the sketch to define a generative model.

  4. Use the generative model to build an estimator.

  5. Profit.

Model “recipes”

  1. Recognize a set of variables to work with. (Data and parameters.)
  2. Define each variable either in terms of the other variables OR in terms of a probability distribution.
  3. The combination of variables and their probability distributions defines a joint generative model that can be used to simulate hypothetical observations and analyze real ones.

Here’s an example:

\[\begin{align*} y_i &\sim \text{Normal}(\mu_i,\sigma) \\ \mu_i &= \beta x_i \\ \beta &\sim \text{Normal}(0,10) \\ \sigma &\sim \text{Exponential}(1) \\ x_i &\sim \text{Normal}(0,1) \\ \end{align*}\]

Model for globe-tossing

Here’s the model for last week’s globe-tossing experiment:

\[\begin{align*} W &\sim \text{Binomial}(N,p) \\ p &\sim \text{Uniform}(0,1) \\ \end{align*}\]

  • \(W\) is the observed count of water.
  • \(N\) is the total number of tosses.
  • \(p\) is the proportion of water on the globe.

The whole model can be read as:

The count \(W\) is distributed binomially with sample size \(N\) and probability \(p\). The prior for \(p\) is assumed to be uniform between 0 and 1.

Model for globe-tossing

Here’s the model for last week’s globe-tossing experiment:

\[\begin{align*} W &\sim \text{Binomial}(N,p) \\ p &\sim \text{Uniform}(0,1) \\ \end{align*}\]

Estimating the posterior using brms

Last week, we used grid approximation to estimate the posterior distribution. In the video you watched for today, McElreath moves on to something called QUADRATIC APPROXIMATION. It’s good to understand what that’s doing, but you and I are moving right along to MCMC. We won’t go into the details of what’s happening for a few weeks, but let’s start with the code to estimate the model.

Let’s say we tossed the globe 9 times and observed 6 waters:

m1 <-
  brm(data = list(w = 6),
      family = binomial(link = "identity"),
      w | trials(9) ~ 0 + Intercept,
      prior(uniform(0, 1), class = b),
      iter = 5000, warmup = 1000, seed = 3, chains=1,
      file = here("files/models/m21.1"))
1
Data can be a data frame or a list. Just make sure the variable names match what’s in your formula.
2
How you assume your outcome variable is distributed.
3
The formula for your outcome.
4
Priors for every parameter in your model. Here, we only have one parameter, so we only need 1 prior. This happens to be a flat prior between 0 and 1.
5
Some choices about how we want our model to run. We’ll go more into this later.
6
These models can take a long time to run. You can automatically save the output to a file; when you do this, the next time you run this code, it won’t actually estimate the model, but will instead pull the output from your stated file. Be WARNED: if you change the data or the model code, it will NOT restimate your model until you delete the file.

Sampling from the posterior

Grid approximation gave us the calculated probability of each possible value of our parameter, \(p\). But our method of conducting bayes will no longer give us such a neat solution. Here’s how you get the posterior distribution for \(p\):

samples_from_post = as_draws_df(m1)
samples_from_post
# A draws_df: 4000 iterations, 1 chains, and 3 variables
   b_Intercept lprior lp__
1         0.46      0 -2.1
2         0.49      0 -1.9
3         0.56      0 -1.5
4         0.66      0 -1.3
5         0.71      0 -1.3
6         0.62      0 -1.3
7         0.68      0 -1.3
8         0.55      0 -1.5
9         0.63      0 -1.3
10        0.67      0 -1.3
# ... with 3990 more draws
# ... hidden reserved variables {'.chain', '.iteration', '.draw'}
samples_from_post %>%  
  ggplot(aes(x=b_Intercept)) +
  geom_density(fill = "grey", color = "white") +
  labs(x="Proportion water")

Posterior predictive distribution

ppd = posterior_predict(m1)
dim(ppd)
[1] 4000    1
ppd
        [,1]
   [1,]    3
   [2,]    2
   [3,]    6
   [4,]    6
   [5,]    5
   [6,]    6
   [7,]    2
   [8,]    4
   [9,]    8
  [10,]    7
  [11,]    5
  [12,]    8
  [13,]    5
  [14,]    8
  [15,]    8
  [16,]    8
  [17,]    8
  [18,]    9
  [19,]    7
  [20,]    6
  [21,]    9
  [22,]    8
  [23,]    3
  [24,]    7
  [25,]    4
  [26,]    8
  [27,]    9
  [28,]    4
  [29,]    6
  [30,]    5
  [31,]    6
  [32,]    2
  [33,]    7
  [34,]    7
  [35,]    5
  [36,]    5
  [37,]    7
  [38,]    3
  [39,]    6
  [40,]    6
  [41,]    5
  [42,]    8
  [43,]    7
  [44,]    5
  [45,]    7
  [46,]    7
  [47,]    9
  [48,]    4
  [49,]    3
  [50,]    8
  [51,]    6
  [52,]    7
  [53,]    5
  [54,]    5
  [55,]    5
  [56,]    5
  [57,]    4
  [58,]    3
  [59,]    4
  [60,]    5
  [61,]    7
  [62,]    8
  [63,]    6
  [64,]    5
  [65,]    4
  [66,]    7
  [67,]    3
  [68,]    8
  [69,]    7
  [70,]    9
  [71,]    6
  [72,]    4
  [73,]    5
  [74,]    7
  [75,]    4
  [76,]    9
  [77,]    6
  [78,]    8
  [79,]    8
  [80,]    8
  [81,]    6
  [82,]    3
  [83,]    6
  [84,]    7
  [85,]    7
  [86,]    5
  [87,]    6
  [88,]    5
  [89,]    7
  [90,]    6
  [91,]    8
  [92,]    5
  [93,]    5
  [94,]    5
  [95,]    5
  [96,]    9
  [97,]    7
  [98,]    7
  [99,]    6
 [100,]    8
 [101,]    4
 [102,]    7
 [103,]    5
 [104,]    6
 [105,]    8
 [106,]    8
 [107,]    8
 [108,]    7
 [109,]    6
 [110,]    5
 [111,]    3
 [112,]    5
 [113,]    4
 [114,]    6
 [115,]    5
 [116,]    6
 [117,]    6
 [118,]    3
 [119,]    7
 [120,]    4
 [121,]    6
 [122,]    6
 [123,]    7
 [124,]    7
 [125,]    9
 [126,]    6
 [127,]    8
 [128,]    4
 [129,]    8
 [130,]    6
 [131,]    9
 [132,]    7
 [133,]    4
 [134,]    9
 [135,]    6
 [136,]    4
 [137,]    7
 [138,]    6
 [139,]    6
 [140,]    4
 [141,]    8
 [142,]    5
 [143,]    6
 [144,]    8
 [145,]    6
 [146,]    7
 [147,]    6
 [148,]    6
 [149,]    4
 [150,]    6
 [151,]    4
 [152,]    6
 [153,]    8
 [154,]    7
 [155,]    7
 [156,]    5
 [157,]    3
 [158,]    8
 [159,]    6
 [160,]    3
 [161,]    8
 [162,]    6
 [163,]    2
 [164,]    5
 [165,]    4
 [166,]    6
 [167,]    4
 [168,]    7
 [169,]    7
 [170,]    8
 [171,]    7
 [172,]    2
 [173,]    5
 [174,]    8
 [175,]    6
 [176,]    4
 [177,]    8
 [178,]    7
 [179,]    8
 [180,]    6
 [181,]    6
 [182,]    7
 [183,]    8
 [184,]    7
 [185,]    7
 [186,]    7
 [187,]    7
 [188,]    8
 [189,]    4
 [190,]    4
 [191,]    6
 [192,]    8
 [193,]    2
 [194,]    7
 [195,]    4
 [196,]    6
 [197,]    5
 [198,]    7
 [199,]    7
 [200,]    8
 [201,]    9
 [202,]    8
 [203,]    8
 [204,]    7
 [205,]    6
 [206,]    8
 [207,]    5
 [208,]    4
 [209,]    7
 [210,]    6
 [211,]    8
 [212,]    7
 [213,]    8
 [214,]    9
 [215,]    8
 [216,]    7
 [217,]    8
 [218,]    7
 [219,]    3
 [220,]    4
 [221,]    6
 [222,]    4
 [223,]    4
 [224,]    7
 [225,]    6
 [226,]    6
 [227,]    4
 [228,]    8
 [229,]    8
 [230,]    8
 [231,]    7
 [232,]    5
 [233,]    6
 [234,]    6
 [235,]    6
 [236,]    6
 [237,]    6
 [238,]    7
 [239,]    6
 [240,]    4
 [241,]    4
 [242,]    2
 [243,]    4
 [244,]    8
 [245,]    8
 [246,]    9
 [247,]    7
 [248,]    6
 [249,]    5
 [250,]    4
 [251,]    3
 [252,]    7
 [253,]    5
 [254,]    4
 [255,]    4
 [256,]    8
 [257,]    4
 [258,]    6
 [259,]    3
 [260,]    4
 [261,]    8
 [262,]    6
 [263,]    6
 [264,]    7
 [265,]    8
 [266,]    6
 [267,]    5
 [268,]    6
 [269,]    6
 [270,]    5
 [271,]    6
 [272,]    3
 [273,]    5
 [274,]    6
 [275,]    4
 [276,]    6
 [277,]    8
 [278,]    9
 [279,]    9
 [280,]    8
 [281,]    5
 [282,]    8
 [283,]    2
 [284,]    4
 [285,]    6
 [286,]    3
 [287,]    2
 [288,]    9
 [289,]    5
 [290,]    6
 [291,]    7
 [292,]    7
 [293,]    8
 [294,]    5
 [295,]    2
 [296,]    6
 [297,]    9
 [298,]    5
 [299,]    7
 [300,]    3
 [301,]    8
 [302,]    8
 [303,]    6
 [304,]    7
 [305,]    8
 [306,]    7
 [307,]    8
 [308,]    8
 [309,]    5
 [310,]    8
 [311,]    5
 [312,]    6
 [313,]    8
 [314,]    7
 [315,]    7
 [316,]    7
 [317,]    9
 [318,]    6
 [319,]    8
 [320,]    6
 [321,]    3
 [322,]    4
 [323,]    4
 [324,]    2
 [325,]    5
 [326,]    6
 [327,]    4
 [328,]    9
 [329,]    4
 [330,]    8
 [331,]    6
 [332,]    6
 [333,]    9
 [334,]    6
 [335,]    8
 [336,]    6
 [337,]    8
 [338,]    4
 [339,]    2
 [340,]    0
 [341,]    4
 [342,]    6
 [343,]    8
 [344,]    6
 [345,]    7
 [346,]    5
 [347,]    3
 [348,]    9
 [349,]    8
 [350,]    6
 [351,]    7
 [352,]    6
 [353,]    6
 [354,]    6
 [355,]    4
 [356,]    7
 [357,]    6
 [358,]    4
 [359,]    1
 [360,]    3
 [361,]    5
 [362,]    4
 [363,]    7
 [364,]    8
 [365,]    6
 [366,]    7
 [367,]    7
 [368,]    3
 [369,]    7
 [370,]    8
 [371,]    1
 [372,]    3
 [373,]    1
 [374,]    4
 [375,]    3
 [376,]    3
 [377,]    4
 [378,]    4
 [379,]    4
 [380,]    5
 [381,]    4
 [382,]    3
 [383,]    7
 [384,]    8
 [385,]    6
 [386,]    9
 [387,]    6
 [388,]    9
 [389,]    7
 [390,]    8
 [391,]    5
 [392,]    3
 [393,]    3
 [394,]    8
 [395,]    7
 [396,]    8
 [397,]    5
 [398,]    8
 [399,]    9
 [400,]    8
 [401,]    9
 [402,]    5
 [403,]    4
 [404,]    7
 [405,]    8
 [406,]    6
 [407,]    8
 [408,]    8
 [409,]    8
 [410,]    9
 [411,]    5
 [412,]    7
 [413,]    7
 [414,]    2
 [415,]    6
 [416,]    6
 [417,]    7
 [418,]    6
 [419,]    6
 [420,]    7
 [421,]    7
 [422,]    6
 [423,]    8
 [424,]    7
 [425,]    8
 [426,]    5
 [427,]    4
 [428,]    9
 [429,]    8
 [430,]    5
 [431,]    4
 [432,]    6
 [433,]    7
 [434,]    7
 [435,]    9
 [436,]    5
 [437,]    5
 [438,]    4
 [439,]    6
 [440,]    6
 [441,]    6
 [442,]    6
 [443,]    7
 [444,]    7
 [445,]    4
 [446,]    4
 [447,]    5
 [448,]    8
 [449,]    4
 [450,]    6
 [451,]    3
 [452,]    7
 [453,]    6
 [454,]    8
 [455,]    8
 [456,]    6
 [457,]    6
 [458,]    5
 [459,]    3
 [460,]    5
 [461,]    4
 [462,]    8
 [463,]    7
 [464,]    6
 [465,]    9
 [466,]    6
 [467,]    5
 [468,]    7
 [469,]    5
 [470,]    3
 [471,]    6
 [472,]    4
 [473,]    8
 [474,]    6
 [475,]    5
 [476,]    3
 [477,]    7
 [478,]    6
 [479,]    5
 [480,]    5
 [481,]    6
 [482,]    3
 [483,]    3
 [484,]    4
 [485,]    5
 [486,]    4
 [487,]    6
 [488,]    3
 [489,]    6
 [490,]    5
 [491,]    8
 [492,]    8
 [493,]    3
 [494,]    4
 [495,]    7
 [496,]    6
 [497,]    7
 [498,]    7
 [499,]    6
 [500,]    2
 [501,]    3
 [502,]    5
 [503,]    8
 [504,]    7
 [505,]    7
 [506,]    5
 [507,]    7
 [508,]    8
 [509,]    6
 [510,]    4
 [511,]    7
 [512,]    5
 [513,]    5
 [514,]    6
 [515,]    7
 [516,]    8
 [517,]    3
 [518,]    6
 [519,]    6
 [520,]    3
 [521,]    3
 [522,]    3
 [523,]    2
 [524,]    5
 [525,]    5
 [526,]    8
 [527,]    8
 [528,]    6
 [529,]    5
 [530,]    5
 [531,]    7
 [532,]    7
 [533,]    5
 [534,]    8
 [535,]    5
 [536,]    2
 [537,]    7
 [538,]    5
 [539,]    8
 [540,]    5
 [541,]    4
 [542,]    6
 [543,]    6
 [544,]    6
 [545,]    3
 [546,]    7
 [547,]    5
 [548,]    4
 [549,]    5
 [550,]    5
 [551,]    7
 [552,]    9
 [553,]    4
 [554,]    6
 [555,]    5
 [556,]    5
 [557,]    3
 [558,]    1
 [559,]    3
 [560,]    7
 [561,]    3
 [562,]    5
 [563,]    7
 [564,]    8
 [565,]    7
 [566,]    2
 [567,]    2
 [568,]    7
 [569,]    7
 [570,]    8
 [571,]    5
 [572,]    5
 [573,]    4
 [574,]    2
 [575,]    6
 [576,]    4
 [577,]    6
 [578,]    8
 [579,]    5
 [580,]    4
 [581,]    5
 [582,]    6
 [583,]    3
 [584,]    6
 [585,]    4
 [586,]    3
 [587,]    2
 [588,]    3
 [589,]    2
 [590,]    3
 [591,]    4
 [592,]    5
 [593,]    7
 [594,]    5
 [595,]    6
 [596,]    7
 [597,]    5
 [598,]    2
 [599,]    6
 [600,]    4
 [601,]    7
 [602,]    5
 [603,]    6
 [604,]    6
 [605,]    9
 [606,]    6
 [607,]    7
 [608,]    4
 [609,]    3
 [610,]    5
 [611,]    4
 [612,]    8
 [613,]    5
 [614,]    8
 [615,]    7
 [616,]    1
 [617,]    3
 [618,]    8
 [619,]    8
 [620,]    3
 [621,]    7
 [622,]    7
 [623,]    9
 [624,]    6
 [625,]    4
 [626,]    4
 [627,]    2
 [628,]    5
 [629,]    6
 [630,]    5
 [631,]    5
 [632,]    4
 [633,]    5
 [634,]    5
 [635,]    3
 [636,]    9
 [637,]    6
 [638,]    4
 [639,]    6
 [640,]    3
 [641,]    6
 [642,]    6
 [643,]    8
 [644,]    7
 [645,]    7
 [646,]    5
 [647,]    5
 [648,]    4
 [649,]    2
 [650,]    5
 [651,]    7
 [652,]    7
 [653,]    6
 [654,]    3
 [655,]    6
 [656,]    7
 [657,]    8
 [658,]    8
 [659,]    7
 [660,]    2
 [661,]    6
 [662,]    2
 [663,]    4
 [664,]    3
 [665,]    8
 [666,]    7
 [667,]    9
 [668,]    5
 [669,]    7
 [670,]    7
 [671,]    5
 [672,]    7
 [673,]    5
 [674,]    3
 [675,]    3
 [676,]    7
 [677,]    7
 [678,]    3
 [679,]    8
 [680,]    7
 [681,]    6
 [682,]    4
 [683,]    3
 [684,]    5
 [685,]    7
 [686,]    9
 [687,]    7
 [688,]    3
 [689,]    6
 [690,]    5
 [691,]    6
 [692,]    5
 [693,]    3
 [694,]    2
 [695,]    4
 [696,]    9
 [697,]    5
 [698,]    2
 [699,]    3
 [700,]    3
 [701,]    8
 [702,]    7
 [703,]    7
 [704,]    5
 [705,]    2
 [706,]    7
 [707,]    4
 [708,]    5
 [709,]    8
 [710,]    8
 [711,]    3
 [712,]    3
 [713,]    6
 [714,]    7
 [715,]    6
 [716,]    6
 [717,]    7
 [718,]    8
 [719,]    4
 [720,]    3
 [721,]    6
 [722,]    6
 [723,]    3
 [724,]    8
 [725,]    7
 [726,]    7
 [727,]    3
 [728,]    6
 [729,]    5
 [730,]    3
 [731,]    5
 [732,]    5
 [733,]    4
 [734,]    7
 [735,]    8
 [736,]    8
 [737,]    6
 [738,]    6
 [739,]    5
 [740,]    7
 [741,]    7
 [742,]    8
 [743,]    6
 [744,]    7
 [745,]    7
 [746,]    5
 [747,]    8
 [748,]    3
 [749,]    4
 [750,]    5
 [751,]    6
 [752,]    6
 [753,]    6
 [754,]    6
 [755,]    8
 [756,]    6
 [757,]    8
 [758,]    7
 [759,]    7
 [760,]    8
 [761,]    5
 [762,]    7
 [763,]    8
 [764,]    5
 [765,]    7
 [766,]    7
 [767,]    7
 [768,]    6
 [769,]    4
 [770,]    8
 [771,]    6
 [772,]    3
 [773,]    5
 [774,]    4
 [775,]    8
 [776,]    5
 [777,]    3
 [778,]    8
 [779,]    5
 [780,]    3
 [781,]    5
 [782,]    7
 [783,]    6
 [784,]    5
 [785,]    3
 [786,]    3
 [787,]    2
 [788,]    8
 [789,]    3
 [790,]    8
 [791,]    9
 [792,]    7
 [793,]    7
 [794,]    8
 [795,]    3
 [796,]    4
 [797,]    8
 [798,]    6
 [799,]    8
 [800,]    9
 [801,]    9
 [802,]    8
 [803,]    6
 [804,]    5
 [805,]    4
 [806,]    8
 [807,]    5
 [808,]    8
 [809,]    6
 [810,]    7
 [811,]    4
 [812,]    8
 [813,]    6
 [814,]    8
 [815,]    7
 [816,]    7
 [817,]    9
 [818,]    7
 [819,]    6
 [820,]    7
 [821,]    2
 [822,]    6
 [823,]    2
 [824,]    6
 [825,]    5
 [826,]    3
 [827,]    4
 [828,]    3
 [829,]    1
 [830,]    6
 [831,]    7
 [832,]    2
 [833,]    4
 [834,]    6
 [835,]    6
 [836,]    7
 [837,]    7
 [838,]    7
 [839,]    3
 [840,]    3
 [841,]    8
 [842,]    7
 [843,]    7
 [844,]    6
 [845,]    8
 [846,]    5
 [847,]    6
 [848,]    6
 [849,]    5
 [850,]    4
 [851,]    7
 [852,]    8
 [853,]    6
 [854,]    4
 [855,]    8
 [856,]    5
 [857,]    6
 [858,]    6
 [859,]    6
 [860,]    6
 [861,]    9
 [862,]    6
 [863,]    6
 [864,]    6
 [865,]    4
 [866,]    9
 [867,]    7
 [868,]    5
 [869,]    7
 [870,]    5
 [871,]    9
 [872,]    7
 [873,]    6
 [874,]    6
 [875,]    1
 [876,]    2
 [877,]    5
 [878,]    2
 [879,]    1
 [880,]    2
 [881,]    2
 [882,]    8
 [883,]    9
 [884,]    5
 [885,]    3
 [886,]    5
 [887,]    6
 [888,]    5
 [889,]    7
 [890,]    8
 [891,]    5
 [892,]    6
 [893,]    6
 [894,]    8
 [895,]    5
 [896,]    6
 [897,]    6
 [898,]    5
 [899,]    9
 [900,]    7
 [901,]    4
 [902,]    8
 [903,]    4
 [904,]    4
 [905,]    7
 [906,]    8
 [907,]    7
 [908,]    7
 [909,]    8
 [910,]    8
 [911,]    9
 [912,]    9
 [913,]    8
 [914,]    9
 [915,]    4
 [916,]    5
 [917,]    4
 [918,]    6
 [919,]    6
 [920,]    5
 [921,]    5
 [922,]    4
 [923,]    5
 [924,]    3
 [925,]    6
 [926,]    6
 [927,]    3
 [928,]    3
 [929,]    5
 [930,]    3
 [931,]    6
 [932,]    3
 [933,]    7
 [934,]    6
 [935,]    2
 [936,]    3
 [937,]    8
 [938,]    7
 [939,]    8
 [940,]    8
 [941,]    5
 [942,]    3
 [943,]    6
 [944,]    9
 [945,]    6
 [946,]    5
 [947,]    4
 [948,]    7
 [949,]    8
 [950,]    4
 [951,]    5
 [952,]    6
 [953,]    5
 [954,]    5
 [955,]    5
 [956,]    5
 [957,]    4
 [958,]    5
 [959,]    8
 [960,]    8
 [961,]    8
 [962,]    5
 [963,]    6
 [964,]    6
 [965,]    1
 [966,]    2
 [967,]    3
 [968,]    6
 [969,]    4
 [970,]    5
 [971,]    8
 [972,]    8
 [973,]    2
 [974,]    3
 [975,]    5
 [976,]    6
 [977,]    4
 [978,]    7
 [979,]    8
 [980,]    9
 [981,]    8
 [982,]    8
 [983,]    3
 [984,]    8
 [985,]    2
 [986,]    8
 [987,]    6
 [988,]    7
 [989,]    3
 [990,]    3
 [991,]    5
 [992,]    7
 [993,]    6
 [994,]    8
 [995,]    7
 [996,]    7
 [997,]    3
 [998,]    4
 [999,]    6
[1000,]    3
[1001,]    3
[1002,]    7
[1003,]    7
[1004,]    3
[1005,]    6
[1006,]    9
[1007,]    5
[1008,]    8
[1009,]    7
[1010,]    8
[1011,]    7
[1012,]    5
[1013,]    3
[1014,]    5
[1015,]    7
[1016,]    4
[1017,]    7
[1018,]    9
[1019,]    2
[1020,]    7
[1021,]    2
[1022,]    3
[1023,]    1
[1024,]    6
[1025,]    4
[1026,]    7
[1027,]    9
[1028,]    4
[1029,]    5
[1030,]    6
[1031,]    4
[1032,]    9
[1033,]    5
[1034,]    4
[1035,]    5
[1036,]    6
[1037,]    4
[1038,]    6
[1039,]    6
[1040,]    6
[1041,]    3
[1042,]    7
[1043,]    2
[1044,]    5
[1045,]    7
[1046,]    4
[1047,]    6
[1048,]    3
[1049,]    3
[1050,]    7
[1051,]    7
[1052,]    5
[1053,]    7
[1054,]    9
[1055,]    6
[1056,]    6
[1057,]    3
[1058,]    4
[1059,]    4
[1060,]    3
[1061,]    5
[1062,]    7
[1063,]    7
[1064,]    9
[1065,]    7
[1066,]    8
[1067,]    7
[1068,]    9
[1069,]    5
[1070,]    7
[1071,]    8
[1072,]    9
[1073,]    7
[1074,]    5
[1075,]    7
[1076,]    9
[1077,]    8
[1078,]    3
[1079,]    3
[1080,]    6
[1081,]    8
[1082,]    6
[1083,]    6
[1084,]    5
[1085,]    6
[1086,]    5
[1087,]    9
[1088,]    7
[1089,]    8
[1090,]    8
[1091,]    7
[1092,]    8
[1093,]    7
[1094,]    5
[1095,]    7
[1096,]    4
[1097,]    4
[1098,]    7
[1099,]    5
[1100,]    5
[1101,]    6
[1102,]    4
[1103,]    5
[1104,]    7
[1105,]    9
[1106,]    8
[1107,]    7
[1108,]    4
[1109,]    6
[1110,]    6
[1111,]    7
[1112,]    9
[1113,]    5
[1114,]    8
[1115,]    5
[1116,]    7
[1117,]    4
[1118,]    6
[1119,]    6
[1120,]    6
[1121,]    9
[1122,]    5
[1123,]    6
[1124,]    7
[1125,]    8
[1126,]    4
[1127,]    4
[1128,]    5
[1129,]    4
[1130,]    5
[1131,]    2
[1132,]    4
[1133,]    6
[1134,]    7
[1135,]    5
[1136,]    6
[1137,]    6
[1138,]    5
[1139,]    5
[1140,]    3
[1141,]    5
[1142,]    5
[1143,]    6
[1144,]    7
[1145,]    5
[1146,]    8
[1147,]    9
[1148,]    8
[1149,]    7
[1150,]    6
[1151,]    6
[1152,]    6
[1153,]    6
[1154,]    8
[1155,]    9
[1156,]    6
[1157,]    6
[1158,]    4
[1159,]    7
[1160,]    8
[1161,]    5
[1162,]    6
[1163,]    2
[1164,]    4
[1165,]    6
[1166,]    3
[1167,]    4
[1168,]    3
[1169,]    5
[1170,]    7
[1171,]    7
[1172,]    6
[1173,]    3
[1174,]    1
[1175,]    3
[1176,]    8
[1177,]    8
[1178,]    7
[1179,]    6
[1180,]    4
[1181,]    3
[1182,]    4
[1183,]    2
[1184,]    7
[1185,]    4
[1186,]    7
[1187,]    4
[1188,]    2
[1189,]    6
[1190,]    8
[1191,]    8
[1192,]    4
[1193,]    7
[1194,]    6
[1195,]    5
[1196,]    3
[1197,]    7
[1198,]    7
[1199,]    8
[1200,]    4
[1201,]    7
[1202,]    6
[1203,]    5
[1204,]    3
[1205,]    3
[1206,]    5
[1207,]    7
[1208,]    6
[1209,]    6
[1210,]    5
[1211,]    7
[1212,]    7
[1213,]    7
[1214,]    6
[1215,]    6
[1216,]    9
[1217,]    8
[1218,]    9
[1219,]    8
[1220,]    9
[1221,]    5
[1222,]    6
[1223,]    5
[1224,]    5
[1225,]    5
[1226,]    7
[1227,]    9
[1228,]    8
[1229,]    5
[1230,]    6
[1231,]    5
[1232,]    5
[1233,]    5
[1234,]    3
[1235,]    2
[1236,]    4
[1237,]    4
[1238,]    3
[1239,]    7
[1240,]    6
[1241,]    9
[1242,]    8
[1243,]    9
[1244,]    7
[1245,]    6
[1246,]    5
[1247,]    8
[1248,]    6
[1249,]    8
[1250,]    5
[1251,]    7
[1252,]    8
[1253,]    8
[1254,]    6
[1255,]    7
[1256,]    8
[1257,]    9
[1258,]    4
[1259,]    9
[1260,]    6
[1261,]    8
[1262,]    8
[1263,]    6
[1264,]    4
[1265,]    4
[1266,]    4
[1267,]    5
[1268,]    5
[1269,]    6
[1270,]    5
[1271,]    5
[1272,]    6
[1273,]    6
[1274,]    4
[1275,]    7
[1276,]    6
[1277,]    5
[1278,]    7
[1279,]    4
[1280,]    5
[1281,]    5
[1282,]    4
[1283,]    8
[1284,]    5
[1285,]    5
[1286,]    5
[1287,]    4
[1288,]    4
[1289,]    6
[1290,]    5
[1291,]    4
[1292,]    8
[1293,]    6
[1294,]    4
[1295,]    6
[1296,]    8
[1297,]    7
[1298,]    7
[1299,]    8
[1300,]    4
[1301,]    8
[1302,]    7
[1303,]    8
[1304,]    6
[1305,]    8
[1306,]    8
[1307,]    6
[1308,]    7
[1309,]    6
[1310,]    7
[1311,]    5
[1312,]    7
[1313,]    8
[1314,]    8
[1315,]    6
[1316,]    6
[1317,]    2
[1318,]    2
[1319,]    5
[1320,]    7
[1321,]    5
[1322,]    7
[1323,]    6
[1324,]    4
[1325,]    8
[1326,]    7
[1327,]    7
[1328,]    7
[1329,]    9
[1330,]    4
[1331,]    4
[1332,]    4
[1333,]    7
[1334,]    4
[1335,]    4
[1336,]    8
[1337,]    6
[1338,]    7
[1339,]    6
[1340,]    7
[1341,]    8
[1342,]    7
[1343,]    8
[1344,]    9
[1345,]    8
[1346,]    3
[1347,]    4
[1348,]    9
[1349,]    6
[1350,]    3
[1351,]    8
[1352,]    4
[1353,]    3
[1354,]    8
[1355,]    9
[1356,]    9
[1357,]    4
[1358,]    3
[1359,]    5
[1360,]    4
[1361,]    8
[1362,]    5
[1363,]    4
[1364,]    8
[1365,]    5
[1366,]    5
[1367,]    7
[1368,]    6
[1369,]    5
[1370,]    5
[1371,]    6
[1372,]    3
[1373,]    6
[1374,]    6
[1375,]    3
[1376,]    6
[1377,]    3
[1378,]    4
[1379,]    5
[1380,]    5
[1381,]    5
[1382,]    1
[1383,]    4
[1384,]    5
[1385,]    9
[1386,]    4
[1387,]    6
[1388,]    2
[1389,]    5
[1390,]    8
[1391,]    9
[1392,]    5
[1393,]    4
[1394,]    4
[1395,]    4
[1396,]    4
[1397,]    2
[1398,]    6
[1399,]    2
[1400,]    0
[1401,]    6
[1402,]    4
[1403,]    2
[1404,]    3
[1405,]    6
[1406,]    7
[1407,]    8
[1408,]    7
[1409,]    5
[1410,]    8
[1411,]    7
[1412,]    9
[1413,]    3
[1414,]    5
[1415,]    6
[1416,]    7
[1417,]    7
[1418,]    6
[1419,]    7
[1420,]    8
[1421,]    6
[1422,]    6
[1423,]    2
[1424,]    6
[1425,]    7
[1426,]    7
[1427,]    8
[1428,]    8
[1429,]    5
[1430,]    8
[1431,]    4
[1432,]    3
[1433,]    5
[1434,]    6
[1435,]    1
[1436,]    5
[1437,]    1
[1438,]    7
[1439,]    7
[1440,]    6
[1441,]    3
[1442,]    4
[1443,]    5
[1444,]    5
[1445,]    4
[1446,]    7
[1447,]    8
[1448,]    7
[1449,]    8
[1450,]    5
[1451,]    5
[1452,]    9
[1453,]    9
[1454,]    4
[1455,]    6
[1456,]    7
[1457,]    6
[1458,]    6
[1459,]    3
[1460,]    6
[1461,]    7
[1462,]    7
[1463,]    5
[1464,]    6
[1465,]    4
[1466,]    7
[1467,]    6
[1468,]    4
[1469,]    5
[1470,]    7
[1471,]    4
[1472,]    8
[1473,]    4
[1474,]    9
[1475,]    6
[1476,]    6
[1477,]    7
[1478,]    2
[1479,]    7
[1480,]    5
[1481,]    4
[1482,]    6
[1483,]    4
[1484,]    7
[1485,]    6
[1486,]    8
[1487,]    0
[1488,]    4
[1489,]    7
[1490,]    7
[1491,]    5
[1492,]    5
[1493,]    8
[1494,]    7
[1495,]    4
[1496,]    6
[1497,]    8
[1498,]    8
[1499,]    4
[1500,]    5
[1501,]    7
[1502,]    7
[1503,]    5
[1504,]    6
[1505,]    8
[1506,]    5
[1507,]    7
[1508,]    6
[1509,]    7
[1510,]    2
[1511,]    3
[1512,]    7
[1513,]    5
[1514,]    2
[1515,]    8
[1516,]    7
[1517,]    6
[1518,]    6
[1519,]    4
[1520,]    3
[1521,]    3
[1522,]    1
[1523,]    7
[1524,]    8
[1525,]    9
[1526,]    7
[1527,]    6
[1528,]    6
[1529,]    8
[1530,]    7
[1531,]    5
[1532,]    5
[1533,]    5
[1534,]    8
[1535,]    4
[1536,]    5
[1537,]    6
[1538,]    9
[1539,]    4
[1540,]    4
[1541,]    8
[1542,]    6
[1543,]    7
[1544,]    7
[1545,]    3
[1546,]    8
[1547,]    7
[1548,]    5
[1549,]    4
[1550,]    7
[1551,]    8
[1552,]    7
[1553,]    4
[1554,]    3
[1555,]    6
[1556,]    6
[1557,]    7
[1558,]    6
[1559,]    6
[1560,]    5
[1561,]    3
[1562,]    6
[1563,]    5
[1564,]    4
[1565,]    5
[1566,]    5
[1567,]    6
[1568,]    3
[1569,]    4
[1570,]    6
[1571,]    3
[1572,]    3
[1573,]    9
[1574,]    9
[1575,]    8
[1576,]    8
[1577,]    5
[1578,]    2
[1579,]    1
[1580,]    6
[1581,]    4
[1582,]    7
[1583,]    7
[1584,]    5
[1585,]    6
[1586,]    6
[1587,]    4
[1588,]    7
[1589,]    8
[1590,]    6
[1591,]    4
[1592,]    3
[1593,]    4
[1594,]    7
[1595,]    6
[1596,]    6
[1597,]    7
[1598,]    6
[1599,]    7
[1600,]    7
[1601,]    9
[1602,]    6
[1603,]    5
[1604,]    9
[1605,]    7
[1606,]    8
[1607,]    7
[1608,]    7
[1609,]    4
[1610,]    6
[1611,]    7
[1612,]    4
[1613,]    6
[1614,]    3
[1615,]    2
[1616,]    2
[1617,]    7
[1618,]    6
[1619,]    5
[1620,]    6
[1621,]    6
[1622,]    6
[1623,]    6
[1624,]    6
[1625,]    9
[1626,]    7
[1627,]    4
[1628,]    3
[1629,]    3
[1630,]    5
[1631,]    7
[1632,]    7
[1633,]    6
[1634,]    4
[1635,]    3
[1636,]    3
[1637,]    5
[1638,]    5
[1639,]    3
[1640,]    4
[1641,]    5
[1642,]    2
[1643,]    4
[1644,]    7
[1645,]    6
[1646,]    6
[1647,]    5
[1648,]    7
[1649,]    8
[1650,]    5
[1651,]    8
[1652,]    5
[1653,]    4
[1654,]    6
[1655,]    7
[1656,]    4
[1657,]    5
[1658,]    6
[1659,]    9
[1660,]    3
[1661,]    5
[1662,]    5
[1663,]    7
[1664,]    4
[1665,]    5
[1666,]    4
[1667,]    4
[1668,]    3
[1669,]    7
[1670,]    8
[1671,]    3
[1672,]    7
[1673,]    7
[1674,]    7
[1675,]    6
[1676,]    3
[1677,]    9
[1678,]    5
[1679,]    2
[1680,]    5
[1681,]    4
[1682,]    3
[1683,]    5
[1684,]    4
[1685,]    5
[1686,]    4
[1687,]    1
[1688,]    8
[1689,]    7
[1690,]    5
[1691,]    3
[1692,]    6
[1693,]    4
[1694,]    4
[1695,]    9
[1696,]    6
[1697,]    4
[1698,]    6
[1699,]    8
[1700,]    6
[1701,]    8
[1702,]    5
[1703,]    7
[1704,]    2
[1705,]    4
[1706,]    1
[1707,]    7
[1708,]    7
[1709,]    8
[1710,]    5
[1711,]    5
[1712,]    4
[1713,]    6
[1714,]    5
[1715,]    9
[1716,]    8
[1717,]    4
[1718,]    4
[1719,]    7
[1720,]    7
[1721,]    6
[1722,]    7
[1723,]    4
[1724,]    4
[1725,]    4
[1726,]    6
[1727,]    6
[1728,]    6
[1729,]    7
[1730,]    6
[1731,]    3
[1732,]    6
[1733,]    7
[1734,]    8
[1735,]    5
[1736,]    7
[1737,]    6
[1738,]    6
[1739,]    7
[1740,]    5
[1741,]    6
[1742,]    9
[1743,]    6
[1744,]    5
[1745,]    7
[1746,]    8
[1747,]    7
[1748,]    6
[1749,]    7
[1750,]    7
[1751,]    6
[1752,]    4
[1753,]    7
[1754,]    5
[1755,]    7
[1756,]    7
[1757,]    4
[1758,]    3
[1759,]    2
[1760,]    2
[1761,]    1
[1762,]    8
[1763,]    3
[1764,]    5
[1765,]    5
[1766,]    7
[1767,]    6
[1768,]    5
[1769,]    6
[1770,]    6
[1771,]    4
[1772,]    3
[1773,]    6
[1774,]    2
[1775,]    1
[1776,]    4
[1777,]    3
[1778,]    7
[1779,]    4
[1780,]    7
[1781,]    8
[1782,]    5
[1783,]    9
[1784,]    9
[1785,]    9
[1786,]    9
[1787,]    6
[1788,]    5
[1789,]    4
[1790,]    7
[1791,]    6
[1792,]    7
[1793,]    7
[1794,]    6
[1795,]    5
[1796,]    6
[1797,]    7
[1798,]    5
[1799,]    4
[1800,]    4
[1801,]    1
[1802,]    6
[1803,]    7
[1804,]    5
[1805,]    9
[1806,]    7
[1807,]    5
[1808,]    6
[1809,]    6
[1810,]    7
[1811,]    8
[1812,]    6
[1813,]    5
[1814,]    6
[1815,]    6
[1816,]    4
[1817,]    6
[1818,]    7
[1819,]    4
[1820,]    0
[1821,]    4
[1822,]    4
[1823,]    7
[1824,]    6
[1825,]    7
[1826,]    6
[1827,]    6
[1828,]    7
[1829,]    7
[1830,]    8
[1831,]    7
[1832,]    8
[1833,]    9
[1834,]    3
[1835,]    9
[1836,]    1
[1837,]    7
[1838,]    4
[1839,]    8
[1840,]    6
[1841,]    4
[1842,]    7
[1843,]    5
[1844,]    4
[1845,]    5
[1846,]    5
[1847,]    9
[1848,]    6
[1849,]    4
[1850,]    5
[1851,]    6
[1852,]    5
[1853,]    6
[1854,]    5
[1855,]    7
[1856,]    5
[1857,]    8
[1858,]    5
[1859,]    6
[1860,]    6
[1861,]    4
[1862,]    6
[1863,]    5
[1864,]    4
[1865,]    7
[1866,]    5
[1867,]    6
[1868,]    8
[1869,]    7
[1870,]    7
[1871,]    8
[1872,]    5
[1873,]    3
[1874,]    6
[1875,]    6
[1876,]    9
[1877,]    6
[1878,]    8
[1879,]    5
[1880,]    7
[1881,]    7
[1882,]    8
[1883,]    8
[1884,]    8
[1885,]    7
[1886,]    8
[1887,]    7
[1888,]    5
[1889,]    7
[1890,]    6
[1891,]    3
[1892,]    6
[1893,]    9
[1894,]    8
[1895,]    8
[1896,]    8
[1897,]    4
[1898,]    6
[1899,]    7
[1900,]    6
[1901,]    5
[1902,]    7
[1903,]    8
[1904,]    8
[1905,]    7
[1906,]    9
[1907,]    7
[1908,]    7
[1909,]    8
[1910,]    9
[1911,]    7
[1912,]    6
[1913,]    4
[1914,]    7
[1915,]    3
[1916,]    3
[1917,]    4
[1918,]    6
[1919,]    6
[1920,]    6
[1921,]    6
[1922,]    3
[1923,]    5
[1924,]    6
[1925,]    7
[1926,]    7
[1927,]    5
[1928,]    6
[1929,]    7
[1930,]    5
[1931,]    4
[1932,]    6
[1933,]    8
[1934,]    6
[1935,]    5
[1936,]    6
[1937,]    3
[1938,]    5
[1939,]    3
[1940,]    4
[1941,]    5
[1942,]    3
[1943,]    3
[1944,]    1
[1945,]    6
[1946,]    5
[1947,]    8
[1948,]    5
[1949,]    4
[1950,]    4
[1951,]    6
[1952,]    5
[1953,]    5
[1954,]    6
[1955,]    5
[1956,]    8
[1957,]    6
[1958,]    6
[1959,]    8
[1960,]    9
[1961,]    9
[1962,]    3
[1963,]    8
[1964,]    3
[1965,]    6
[1966,]    5
[1967,]    5
[1968,]    7
[1969,]    6
[1970,]    6
[1971,]    6
[1972,]    7
[1973,]    4
[1974,]    3
[1975,]    5
[1976,]    7
[1977,]    6
[1978,]    6
[1979,]    6
[1980,]    5
[1981,]    5
[1982,]    3
[1983,]    5
[1984,]    7
[1985,]    5
[1986,]    7
[1987,]    7
[1988,]    7
[1989,]    6
[1990,]    5
[1991,]    5
[1992,]    3
[1993,]    4
[1994,]    3
[1995,]    4
[1996,]    2
[1997,]    3
[1998,]    2
[1999,]    5
[2000,]    8
[2001,]    6
[2002,]    5
[2003,]    2
[2004,]    3
[2005,]    4
[2006,]    6
[2007,]    7
[2008,]    6
[2009,]    7
[2010,]    4
[2011,]    6
[2012,]    4
[2013,]    9
[2014,]    8
[2015,]    8
[2016,]    8
[2017,]    5
[2018,]    5
[2019,]    5
[2020,]    7
[2021,]    7
[2022,]    5
[2023,]    3
[2024,]    7
[2025,]    7
[2026,]    6
[2027,]    6
[2028,]    4
[2029,]    8
[2030,]    5
[2031,]    9
[2032,]    8
[2033,]    7
[2034,]    2
[2035,]    6
[2036,]    4
[2037,]    7
[2038,]    5
[2039,]    5
[2040,]    5
[2041,]    7
[2042,]    3
[2043,]    7
[2044,]    7
[2045,]    4
[2046,]    7
[2047,]    6
[2048,]    6
[2049,]    2
[2050,]    4
[2051,]    6
[2052,]    5
[2053,]    5
[2054,]    2
[2055,]    3
[2056,]    5
[2057,]    3
[2058,]    2
[2059,]    6
[2060,]    2
[2061,]    9
[2062,]    4
[2063,]    7
[2064,]    8
[2065,]    7
[2066,]    8
[2067,]    5
[2068,]    8
[2069,]    8
[2070,]    6
[2071,]    5
[2072,]    7
[2073,]    6
[2074,]    6
[2075,]    5
[2076,]    5
[2077,]    2
[2078,]    3
[2079,]    2
[2080,]    4
[2081,]    5
[2082,]    5
[2083,]    7
[2084,]    8
[2085,]    6
[2086,]    6
[2087,]    5
[2088,]    5
[2089,]    8
[2090,]    4
[2091,]    8
[2092,]    7
[2093,]    5
[2094,]    2
[2095,]    4
[2096,]    6
[2097,]    8
[2098,]    7
[2099,]    7
[2100,]    8
[2101,]    8
[2102,]    5
[2103,]    3
[2104,]    4
[2105,]    6
[2106,]    8
[2107,]    7
[2108,]    8
[2109,]    9
[2110,]    6
[2111,]    6
[2112,]    8
[2113,]    8
[2114,]    7
[2115,]    7
[2116,]    9
[2117,]    8
[2118,]    6
[2119,]    5
[2120,]    5
[2121,]    5
[2122,]    1
[2123,]    5
[2124,]    3
[2125,]    5
[2126,]    7
[2127,]    8
[2128,]    6
[2129,]    4
[2130,]    2
[2131,]    6
[2132,]    4
[2133,]    4
[2134,]    7
[2135,]    7
[2136,]    5
[2137,]    7
[2138,]    7
[2139,]    5
[2140,]    5
[2141,]    5
[2142,]    6
[2143,]    5
[2144,]    5
[2145,]    8
[2146,]    8
[2147,]    5
[2148,]    8
[2149,]    5
[2150,]    7
[2151,]    8
[2152,]    7
[2153,]    6
[2154,]    7
[2155,]    5
[2156,]    6
[2157,]    5
[2158,]    7
[2159,]    6
[2160,]    7
[2161,]    6
[2162,]    6
[2163,]    4
[2164,]    8
[2165,]    5
[2166,]    2
[2167,]    4
[2168,]    6
[2169,]    9
[2170,]    6
[2171,]    7
[2172,]    2
[2173,]    6
[2174,]    6
[2175,]    4
[2176,]    7
[2177,]    3
[2178,]    7
[2179,]    8
[2180,]    9
[2181,]    5
[2182,]    6
[2183,]    2
[2184,]    2
[2185,]    5
[2186,]    6
[2187,]    6
[2188,]    6
[2189,]    9
[2190,]    9
[2191,]    4
[2192,]    7
[2193,]    4
[2194,]    7
[2195,]    7
[2196,]    6
[2197,]    5
[2198,]    9
[2199,]    8
[2200,]    7
[2201,]    7
[2202,]    5
[2203,]    6
[2204,]    7
[2205,]    6
[2206,]    5
[2207,]    8
[2208,]    3
[2209,]    4
[2210,]    5
[2211,]    8
[2212,]    6
[2213,]    5
[2214,]    6
[2215,]    3
[2216,]    5
[2217,]    4
[2218,]    7
[2219,]    5
[2220,]    5
[2221,]    4
[2222,]    4
[2223,]    3
[2224,]    7
[2225,]    5
[2226,]    4
[2227,]    5
[2228,]    6
[2229,]    3
[2230,]    7
[2231,]    5
[2232,]    4
[2233,]    3
[2234,]    8
[2235,]    6
[2236,]    6
[2237,]    9
[2238,]    8
[2239,]    5
[2240,]    5
[2241,]    9
[2242,]    6
[2243,]    6
[2244,]    7
[2245,]    2
[2246,]    5
[2247,]    7
[2248,]    7
[2249,]    8
[2250,]    8
[2251,]    8
[2252,]    5
[2253,]    3
[2254,]    7
[2255,]    3
[2256,]    4
[2257,]    4
[2258,]    5
[2259,]    6
[2260,]    3
[2261,]    8
[2262,]    7
[2263,]    8
[2264,]    6
[2265,]    7
[2266,]    7
[2267,]    5
[2268,]    6
[2269,]    6
[2270,]    2
[2271,]    7
[2272,]    4
[2273,]    7
[2274,]    7
[2275,]    6
[2276,]    7
[2277,]    8
[2278,]    8
[2279,]    7
[2280,]    7
[2281,]    6
[2282,]    7
[2283,]    5
[2284,]    6
[2285,]    5
[2286,]    5
[2287,]    7
[2288,]    7
[2289,]    6
[2290,]    6
[2291,]    4
[2292,]    7
[2293,]    7
[2294,]    7
[2295,]    8
[2296,]    6
[2297,]    5
[2298,]    6
[2299,]    8
[2300,]    7
[2301,]    8
[2302,]    4
[2303,]    5
[2304,]    7
[2305,]    3
[2306,]    6
[2307,]    7
[2308,]    4
[2309,]    4
[2310,]    5
[2311,]    3
[2312,]    6
[2313,]    5
[2314,]    7
[2315,]    4
[2316,]    5
[2317,]    2
[2318,]    3
[2319,]    4
[2320,]    5
[2321,]    6
[2322,]    8
[2323,]    6
[2324,]    6
[2325,]    8
[2326,]    5
[2327,]    4
[2328,]    9
[2329,]    9
[2330,]    4
[2331,]    9
[2332,]    8
[2333,]    7
[2334,]    8
[2335,]    6
[2336,]    6
[2337,]    5
[2338,]    7
[2339,]    6
[2340,]    5
[2341,]    6
[2342,]    9
[2343,]    8
[2344,]    9
[2345,]    4
[2346,]    7
[2347,]    8
[2348,]    7
[2349,]    2
[2350,]    6
[2351,]    5
[2352,]    7
[2353,]    3
[2354,]    2
[2355,]    1
[2356,]    7
[2357,]    6
[2358,]    6
[2359,]    4
[2360,]    6
[2361,]    7
[2362,]    7
[2363,]    6
[2364,]    6
[2365,]    8
[2366,]    2
[2367,]    3
[2368,]    6
[2369,]    6
[2370,]    6
[2371,]    5
[2372,]    7
[2373,]    5
[2374,]    3
[2375,]    4
[2376,]    7
[2377,]    8
[2378,]    5
[2379,]    4
[2380,]    6
[2381,]    8
[2382,]    7
[2383,]    8
[2384,]    7
[2385,]    4
[2386,]    4
[2387,]    3
[2388,]    4
[2389,]    0
[2390,]    4
[2391,]    7
[2392,]    5
[2393,]    4
[2394,]    6
[2395,]    7
[2396,]    6
[2397,]    3
[2398,]    5
[2399,]    8
[2400,]    7
[2401,]    7
[2402,]    9
[2403,]    4
[2404,]    6
[2405,]    6
[2406,]    5
[2407,]    7
[2408,]    9
[2409,]    8
[2410,]    7
[2411,]    8
[2412,]    7
[2413,]    5
[2414,]    6
[2415,]    6
[2416,]    7
[2417,]    7
[2418,]    7
[2419,]    8
[2420,]    7
[2421,]    4
[2422,]    5
[2423,]    4
[2424,]    5
[2425,]    3
[2426,]    4
[2427,]    4
[2428,]    5
[2429,]    6
[2430,]    7
[2431,]    6
[2432,]    1
[2433,]    4
[2434,]    6
[2435,]    5
[2436,]    9
[2437,]    7
[2438,]    5
[2439,]    4
[2440,]    5
[2441,]    5
[2442,]    4
[2443,]    3
[2444,]    1
[2445,]    7
[2446,]    3
[2447,]    3
[2448,]    5
[2449,]    6
[2450,]    4
[2451,]    7
[2452,]    8
[2453,]    4
[2454,]    7
[2455,]    9
[2456,]    7
[2457,]    8
[2458,]    4
[2459,]    7
[2460,]    7
[2461,]    7
[2462,]    8
[2463,]    6
[2464,]    8
[2465,]    6
[2466,]    8
[2467,]    7
[2468,]    9
[2469,]    6
[2470,]    5
[2471,]    2
[2472,]    8
[2473,]    4
[2474,]    6
[2475,]    7
[2476,]    6
[2477,]    6
[2478,]    6
[2479,]    7
[2480,]    8
[2481,]    7
[2482,]    9
[2483,]    5
[2484,]    6
[2485,]    6
[2486,]    5
[2487,]    7
[2488,]    8
[2489,]    8
[2490,]    4
[2491,]    7
[2492,]    4
[2493,]    5
[2494,]    6
[2495,]    3
[2496,]    5
[2497,]    5
[2498,]    8
[2499,]    6
[2500,]    5
[2501,]    7
[2502,]    6
[2503,]    7
[2504,]    5
[2505,]    8
[2506,]    7
[2507,]    8
[2508,]    7
[2509,]    3
[2510,]    9
[2511,]    4
[2512,]    3
[2513,]    6
[2514,]    8
[2515,]    8
[2516,]    9
[2517,]    4
[2518,]    6
[2519,]    4
[2520,]    5
[2521,]    7
[2522,]    5
[2523,]    8
[2524,]    6
[2525,]    6
[2526,]    5
[2527,]    7
[2528,]    8
[2529,]    7
[2530,]    6
[2531,]    9
[2532,]    8
[2533,]    8
[2534,]    7
[2535,]    5
[2536,]    8
[2537,]    8
[2538,]    6
[2539,]    4
[2540,]    8
[2541,]    6
[2542,]    6
[2543,]    4
[2544,]    3
[2545,]    5
[2546,]    4
[2547,]    6
[2548,]    7
[2549,]    7
[2550,]    2
[2551,]    6
[2552,]    6
[2553,]    6
[2554,]    5
[2555,]    7
[2556,]    6
[2557,]    7
[2558,]    4
[2559,]    7
[2560,]    8
[2561,]    7
[2562,]    6
[2563,]    8
[2564,]    9
[2565,]    9
[2566,]    7
[2567,]    7
[2568,]    5
[2569,]    3
[2570,]    3
[2571,]    7
[2572,]    7
[2573,]    5
[2574,]    7
[2575,]    5
[2576,]    7
[2577,]    5
[2578,]    7
[2579,]    7
[2580,]    5
[2581,]    6
[2582,]    8
[2583,]    9
[2584,]    7
[2585,]    7
[2586,]    6
[2587,]    7
[2588,]    7
[2589,]    7
[2590,]    8
[2591,]    7
[2592,]    5
[2593,]    5
[2594,]    5
[2595,]    4
[2596,]    5
[2597,]    6
[2598,]    6
[2599,]    7
[2600,]    5
[2601,]    7
[2602,]    5
[2603,]    3
[2604,]    6
[2605,]    6
[2606,]    6
[2607,]    6
[2608,]    7
[2609,]    7
[2610,]    6
[2611,]    6
[2612,]    6
[2613,]    5
[2614,]    6
[2615,]    5
[2616,]    4
[2617,]    5
[2618,]    7
[2619,]    6
[2620,]    8
[2621,]    5
[2622,]    5
[2623,]    6
[2624,]    3
[2625,]    7
[2626,]    4
[2627,]    6
[2628,]    7
[2629,]    7
[2630,]    8
[2631,]    7
[2632,]    3
[2633,]    7
[2634,]    7
[2635,]    5
[2636,]    3
[2637,]    6
[2638,]    7
[2639,]    4
[2640,]    9
[2641,]    9
[2642,]    7
[2643,]    7
[2644,]    8
[2645,]    4
[2646,]    4
[2647,]    4
[2648,]    8
[2649,]    9
[2650,]    2
[2651,]    4
[2652,]    2
[2653,]    4
[2654,]    4
[2655,]    3
[2656,]    3
[2657,]    6
[2658,]    6
[2659,]    3
[2660,]    6
[2661,]    6
[2662,]    7
[2663,]    6
[2664,]    7
[2665,]    6
[2666,]    6
[2667,]    8
[2668,]    7
[2669,]    8
[2670,]    9
[2671,]    6
[2672,]    9
[2673,]    3
[2674,]    4
[2675,]    4
[2676,]    4
[2677,]    5
[2678,]    5
[2679,]    9
[2680,]    5
[2681,]    9
[2682,]    7
[2683,]    9
[2684,]    8
[2685,]    6
[2686,]    6
[2687,]    6
[2688,]    4
[2689,]    8
[2690,]    8
[2691,]    6
[2692,]    8
[2693,]    4
[2694,]    8
[2695,]    7
[2696,]    6
[2697,]    2
[2698,]    4
[2699,]    5
[2700,]    7
[2701,]    7
[2702,]    7
[2703,]    4
[2704,]    5
[2705,]    4
[2706,]    5
[2707,]    5
[2708,]    7
[2709,]    6
[2710,]    8
[2711,]    2
[2712,]    3
[2713,]    8
[2714,]    4
[2715,]    9
[2716,]    7
[2717,]    7
[2718,]    8
[2719,]    3
[2720,]    5
[2721,]    7
[2722,]    6
[2723,]    5
[2724,]    4
[2725,]    4
[2726,]    9
[2727,]    9
[2728,]    8
[2729,]    8
[2730,]    5
[2731,]    2
[2732,]    7
[2733,]    5
[2734,]    7
[2735,]    8
[2736,]    6
[2737,]    6
[2738,]    6
[2739,]    7
[2740,]    5
[2741,]    6
[2742,]    4
[2743,]    8
[2744,]    6
[2745,]    4
[2746,]    8
[2747,]    9
[2748,]    8
[2749,]    9
[2750,]    5
[2751,]    8
[2752,]    5
[2753,]    5
[2754,]    6
[2755,]    6
[2756,]    9
[2757,]    8
[2758,]    6
[2759,]    9
[2760,]    9
[2761,]    7
[2762,]    7
[2763,]    9
[2764,]    8
[2765,]    6
[2766,]    6
[2767,]    4
[2768,]    4
[2769,]    3
[2770,]    4
[2771,]    4
[2772,]    6
[2773,]    6
[2774,]    4
[2775,]    4
[2776,]    6
[2777,]    7
[2778,]    7
[2779,]    7
[2780,]    6
[2781,]    3
[2782,]    2
[2783,]    6
[2784,]    4
[2785,]    4
[2786,]    7
[2787,]    3
[2788,]    5
[2789,]    4
[2790,]    7
[2791,]    5
[2792,]    7
[2793,]    6
[2794,]    3
[2795,]    2
[2796,]    7
[2797,]    6
[2798,]    8
[2799,]    8
[2800,]    8
[2801,]    8
[2802,]    6
[2803,]    3
[2804,]    3
[2805,]    5
[2806,]    4
[2807,]    6
[2808,]    5
[2809,]    5
[2810,]    2
[2811,]    4
[2812,]    6
[2813,]    5
[2814,]    9
[2815,]    6
[2816,]    8
[2817,]    9
[2818,]    8
[2819,]    5
[2820,]    9
[2821,]    4
[2822,]    9
[2823,]    8
[2824,]    7
[2825,]    6
[2826,]    5
[2827,]    2
[2828,]    5
[2829,]    2
[2830,]    5
[2831,]    4
[2832,]    7
[2833,]    8
[2834,]    6
[2835,]    7
[2836,]    8
[2837,]    8
[2838,]    7
[2839,]    6
[2840,]    4
[2841,]    5
[2842,]    5
[2843,]    5
[2844,]    6
[2845,]    4
[2846,]    4
[2847,]    4
[2848,]    5
[2849,]    6
[2850,]    3
[2851,]    7
[2852,]    8
[2853,]    8
[2854,]    6
[2855,]    5
[2856,]    4
[2857,]    4
[2858,]    4
[2859,]    5
[2860,]    6
[2861,]    2
[2862,]    5
[2863,]    5
[2864,]    4
[2865,]    4
[2866,]    5
[2867,]    7
[2868,]    8
[2869,]    5
[2870,]    5
[2871,]    4
[2872,]    6
[2873,]    3
[2874,]    3
[2875,]    4
[2876,]    3
[2877,]    5
[2878,]    3
[2879,]    5
[2880,]    8
[2881,]    6
[2882,]    2
[2883,]    8
[2884,]    9
[2885,]    6
[2886,]    7
[2887,]    6
[2888,]    8
[2889,]    5
[2890,]    7
[2891,]    7
[2892,]    5
[2893,]    4
[2894,]    6
[2895,]    5
[2896,]    5
[2897,]    7
[2898,]    4
[2899,]    7
[2900,]    6
[2901,]    7
[2902,]    8
[2903,]    7
[2904,]    5
[2905,]    9
[2906,]    7
[2907,]    6
[2908,]    7
[2909,]    1
[2910,]    1
[2911,]    3
[2912,]    3
[2913,]    2
[2914,]    4
[2915,]    3
[2916,]    8
[2917,]    9
[2918,]    4
[2919,]    7
[2920,]    8
[2921,]    9
[2922,]    6
[2923,]    6
[2924,]    7
[2925,]    6
[2926,]    3
[2927,]    7
[2928,]    3
[2929,]    6
[2930,]    7
[2931,]    9
[2932,]    7
[2933,]    6
[2934,]    5
[2935,]    7
[2936,]    6
[2937,]    6
[2938,]    5
[2939,]    6
[2940,]    4
[2941,]    9
[2942,]    4
[2943,]    4
[2944,]    3
[2945,]    7
[2946,]    4
[2947,]    1
[2948,]    3
[2949,]    8
[2950,]    8
[2951,]    2
[2952,]    4
[2953,]    6
[2954,]    7
[2955,]    8
[2956,]    5
[2957,]    6
[2958,]    3
[2959,]    5
[2960,]    6
[2961,]    4
[2962,]    6
[2963,]    4
[2964,]    8
[2965,]    9
[2966,]    7
[2967,]    8
[2968,]    8
[2969,]    7
[2970,]    8
[2971,]    7
[2972,]    7
[2973,]    4
[2974,]    9
[2975,]    7
[2976,]    8
[2977,]    5
[2978,]    5
[2979,]    6
[2980,]    5
[2981,]    4
[2982,]    6
[2983,]    5
[2984,]    6
[2985,]    7
[2986,]    6
[2987,]    5
[2988,]    7
[2989,]    7
[2990,]    3
[2991,]    6
[2992,]    7
[2993,]    5
[2994,]    7
[2995,]    4
[2996,]    3
[2997,]    3
[2998,]    7
[2999,]    5
[3000,]    4
[3001,]    4
[3002,]    3
[3003,]    9
[3004,]    8
[3005,]    8
[3006,]    4
[3007,]    5
[3008,]    7
[3009,]    1
[3010,]    4
[3011,]    5
[3012,]    4
[3013,]    6
[3014,]    9
[3015,]    8
[3016,]    6
[3017,]    5
[3018,]    8
[3019,]    5
[3020,]    5
[3021,]    8
[3022,]    8
[3023,]    5
[3024,]    6
[3025,]    4
[3026,]    3
[3027,]    6
[3028,]    8
[3029,]    6
[3030,]    7
[3031,]    7
[3032,]    6
[3033,]    5
[3034,]    6
[3035,]    6
[3036,]    3
[3037,]    6
[3038,]    7
[3039,]    7
[3040,]    5
[3041,]    5
[3042,]    4
[3043,]    3
[3044,]    7
[3045,]    4
[3046,]    8
[3047,]    7
[3048,]    6
[3049,]    5
[3050,]    8
[3051,]    3
[3052,]    6
[3053,]    8
[3054,]    4
[3055,]    8
[3056,]    6
[3057,]    6
[3058,]    5
[3059,]    3
[3060,]    1
[3061,]    3
[3062,]    6
[3063,]    8
[3064,]    3
[3065,]    8
[3066,]    8
[3067,]    7
[3068,]    5
[3069,]    6
[3070,]    5
[3071,]    6
[3072,]    6
[3073,]    8
[3074,]    9
[3075,]    7
[3076,]    7
[3077,]    7
[3078,]    4
[3079,]    7
[3080,]    4
[3081,]    6
[3082,]    3
[3083,]    7
[3084,]    5
[3085,]    3
[3086,]    5
[3087,]    6
[3088,]    8
[3089,]    7
[3090,]    8
[3091,]    3
[3092,]    8
[3093,]    3
[3094,]    6
[3095,]    8
[3096,]    7
[3097,]    8
[3098,]    7
[3099,]    8
[3100,]    4
[3101,]    5
[3102,]    2
[3103,]    8
[3104,]    8
[3105,]    8
[3106,]    5
[3107,]    7
[3108,]    6
[3109,]    9
[3110,]    5
[3111,]    7
[3112,]    7
[3113,]    7
[3114,]    9
[3115,]    5
[3116,]    4
[3117,]    6
[3118,]    5
[3119,]    6
[3120,]    5
[3121,]    5
[3122,]    8
[3123,]    3
[3124,]    0
[3125,]    2
[3126,]    4
[3127,]    6
[3128,]    4
[3129,]    2
[3130,]    2
[3131,]    2
[3132,]    4
[3133,]    4
[3134,]    0
[3135,]    3
[3136,]    6
[3137,]    4
[3138,]    4
[3139,]    2
[3140,]    6
[3141,]    6
[3142,]    6
[3143,]    7
[3144,]    2
[3145,]    6
[3146,]    3
[3147,]    6
[3148,]    5
[3149,]    7
[3150,]    9
[3151,]    6
[3152,]    7
[3153,]    6
[3154,]    6
[3155,]    5
[3156,]    7
[3157,]    3
[3158,]    2
[3159,]    3
[3160,]    5
[3161,]    7
[3162,]    5
[3163,]    6
[3164,]    4
[3165,]    3
[3166,]    7
[3167,]    6
[3168,]    5
[3169,]    7
[3170,]    6
[3171,]    6
[3172,]    6
[3173,]    8
[3174,]    7
[3175,]    1
[3176,]    5
[3177,]    6
[3178,]    4
[3179,]    3
[3180,]    8
[3181,]    5
[3182,]    9
[3183,]    5
[3184,]    7
[3185,]    7
[3186,]    8
[3187,]    2
[3188,]    8
[3189,]    7
[3190,]    5
[3191,]    7
[3192,]    8
[3193,]    6
[3194,]    2
[3195,]    5
[3196,]    7
[3197,]    5
[3198,]    6
[3199,]    6
[3200,]    6
[3201,]    5
[3202,]    5
[3203,]    8
[3204,]    6
[3205,]    7
[3206,]    4
[3207,]    3
[3208,]    5
[3209,]    3
[3210,]    2
[3211,]    5
[3212,]    4
[3213,]    3
[3214,]    3
[3215,]    9
[3216,]    7
[3217,]    8
[3218,]    6
[3219,]    5
[3220,]    6
[3221,]    2
[3222,]    6
[3223,]    5
[3224,]    5
[3225,]    6
[3226,]    5
[3227,]    7
[3228,]    3
[3229,]    6
[3230,]    7
[3231,]    1
[3232,]    5
[3233,]    8
[3234,]    7
[3235,]    4
[3236,]    5
[3237,]    4
[3238,]    6
[3239,]    7
[3240,]    7
[3241,]    5
[3242,]    6
[3243,]    3
[3244,]    5
[3245,]    8
[3246,]    2
[3247,]    8
[3248,]    6
[3249,]    5
[3250,]    5
[3251,]    4
[3252,]    9
[3253,]    8
[3254,]    7
[3255,]    5
[3256,]    5
[3257,]    4
[3258,]    4
[3259,]    4
[3260,]    7
[3261,]    6
[3262,]    6
[3263,]    6
[3264,]    5
[3265,]    7
[3266,]    5
[3267,]    7
[3268,]    7
[3269,]    1
[3270,]    4
[3271,]    1
[3272,]    4
[3273,]    4
[3274,]    6
[3275,]    8
[3276,]    6
[3277,]    5
[3278,]    3
[3279,]    3
[3280,]    3
[3281,]    4
[3282,]    7
[3283,]    4
[3284,]    8
[3285,]    8
[3286,]    6
[3287,]    3
[3288,]    3
[3289,]    3
[3290,]    2
[3291,]    8
[3292,]    8
[3293,]    5
[3294,]    6
[3295,]    3
[3296,]    8
[3297,]    5
[3298,]    6
[3299,]    7
[3300,]    7
[3301,]    4
[3302,]    7
[3303,]    5
[3304,]    5
[3305,]    6
[3306,]    7
[3307,]    5
[3308,]    9
[3309,]    6
[3310,]    7
[3311,]    6
[3312,]    4
[3313,]    8
[3314,]    8
[3315,]    5
[3316,]    7
[3317,]    6
[3318,]    5
[3319,]    9
[3320,]    7
[3321,]    8
[3322,]    3
[3323,]    6
[3324,]    6
[3325,]    5
[3326,]    4
[3327,]    7
[3328,]    6
[3329,]    4
[3330,]    7
[3331,]    6
[3332,]    6
[3333,]    6
[3334,]    6
[3335,]    8
[3336,]    8
[3337,]    4
[3338,]    7
[3339,]    7
[3340,]    7
[3341,]    7
[3342,]    3
[3343,]    4
[3344,]    6
[3345,]    8
[3346,]    5
[3347,]    6
[3348,]    7
[3349,]    9
[3350,]    5
[3351,]    4
[3352,]    6
[3353,]    4
[3354,]    5
[3355,]    6
[3356,]    5
[3357,]    3
[3358,]    7
[3359,]    7
[3360,]    8
[3361,]    5
[3362,]    6
[3363,]    7
[3364,]    5
[3365,]    3
[3366,]    4
[3367,]    5
[3368,]    8
[3369,]    7
[3370,]    5
[3371,]    8
[3372,]    8
[3373,]    8
[3374,]    8
[3375,]    9
[3376,]    8
[3377,]    8
[3378,]    7
[3379,]    3
[3380,]    6
[3381,]    6
[3382,]    5
[3383,]    6
[3384,]    6
[3385,]    5
[3386,]    5
[3387,]    6
[3388,]    5
[3389,]    3
[3390,]    6
[3391,]    9
[3392,]    2
[3393,]    4
[3394,]    9
[3395,]    8
[3396,]    6
[3397,]    9
[3398,]    8
[3399,]    8
[3400,]    5
[3401,]    6
[3402,]    6
[3403,]    8
[3404,]    7
[3405,]    5
[3406,]    5
[3407,]    5
[3408,]    7
[3409,]    7
[3410,]    5
[3411,]    5
[3412,]    4
[3413,]    5
[3414,]    4
[3415,]    6
[3416,]    5
[3417,]    5
[3418,]    8
[3419,]    5
[3420,]    8
[3421,]    8
[3422,]    7
[3423,]    5
[3424,]    6
[3425,]    9
[3426,]    9
[3427,]    8
[3428,]    9
[3429,]    8
[3430,]    7
[3431,]    8
[3432,]    7
[3433,]    6
[3434,]    7
[3435,]    8
[3436,]    6
[3437,]    7
[3438,]    4
[3439,]    3
[3440,]    6
[3441,]    4
[3442,]    4
[3443,]    7
[3444,]    7
[3445,]    7
[3446,]    8
[3447,]    5
[3448,]    5
[3449,]    4
[3450,]    8
[3451,]    7
[3452,]    6
[3453,]    8
[3454,]    7
[3455,]    5
[3456,]    6
[3457,]    4
[3458,]    4
[3459,]    7
[3460,]    7
[3461,]    6
[3462,]    7
[3463,]    7
[3464,]    3
[3465,]    2
[3466,]    3
[3467,]    0
[3468,]    4
[3469,]    6
[3470,]    2
[3471,]    4
[3472,]    4
[3473,]    5
[3474,]    5
[3475,]    7
[3476,]    6
[3477,]    6
[3478,]    8
[3479,]    7
[3480,]    4
[3481,]    5
[3482,]    4
[3483,]    8
[3484,]    9
[3485,]    5
[3486,]    5
[3487,]    5
[3488,]    6
[3489,]    7
[3490,]    4
[3491,]    7
[3492,]    7
[3493,]    5
[3494,]    5
[3495,]    4
[3496,]    2
[3497,]    4
[3498,]    6
[3499,]    6
[3500,]    8
[3501,]    5
[3502,]    7
[3503,]    6
[3504,]    6
[3505,]    9
[3506,]    6
[3507,]    7
[3508,]    8
[3509,]    8
[3510,]    3
[3511,]    7
[3512,]    8
[3513,]    4
[3514,]    6
[3515,]    5
[3516,]    7
[3517,]    3
[3518,]    4
[3519,]    8
[3520,]    5
[3521,]    7
[3522,]    7
[3523,]    9
[3524,]    7
[3525,]    8
[3526,]    7
[3527,]    3
[3528,]    3
[3529,]    4
[3530,]    4
[3531,]    3
[3532,]    1
[3533,]    4
[3534,]    6
[3535,]    5
[3536,]    5
[3537,]    7
[3538,]    6
[3539,]    5
[3540,]    6
[3541,]    6
[3542,]    4
[3543,]    9
[3544,]    3
[3545,]    4
[3546,]    2
[3547,]    1
[3548,]    6
[3549,]    8
[3550,]    7
[3551,]    9
[3552,]    7
[3553,]    7
[3554,]    4
[3555,]    6
[3556,]    3
[3557,]    7
[3558,]    6
[3559,]    9
[3560,]    7
[3561,]    9
[3562,]    7
[3563,]    7
[3564,]    8
[3565,]    8
[3566,]    7
[3567,]    5
[3568,]    3
[3569,]    7
[3570,]    7
[3571,]    4
[3572,]    6
[3573,]    4
[3574,]    6
[3575,]    6
[3576,]    6
[3577,]    6
[3578,]    7
[3579,]    6
[3580,]    8
[3581,]    8
[3582,]    3
[3583,]    4
[3584,]    5
[3585,]    6
[3586,]    9
[3587,]    7
[3588,]    4
[3589,]    5
[3590,]    3
[3591,]    6
[3592,]    7
[3593,]    8
[3594,]    7
[3595,]    7
[3596,]    5
[3597,]    5
[3598,]    6
[3599,]    5
[3600,]    1
[3601,]    9
[3602,]    6
[3603,]    8
[3604,]    7
[3605,]    8
[3606,]    8
[3607,]    5
[3608,]    5
[3609,]    7
[3610,]    5
[3611,]    2
[3612,]    5
[3613,]    6
[3614,]    4
[3615,]    6
[3616,]    2
[3617,]    8
[3618,]    4
[3619,]    8
[3620,]    5
[3621,]    6
[3622,]    5
[3623,]    5
[3624,]    6
[3625,]    6
[3626,]    7
[3627,]    7
[3628,]    6
[3629,]    7
[3630,]    5
[3631,]    5
[3632,]    5
[3633,]    6
[3634,]    8
[3635,]    3
[3636,]    6
[3637,]    7
[3638,]    7
[3639,]    4
[3640,]    7
[3641,]    4
[3642,]    4
[3643,]    6
[3644,]    6
[3645,]    6
[3646,]    9
[3647,]    7
[3648,]    6
[3649,]    7
[3650,]    4
[3651,]    7
[3652,]    8
[3653,]    5
[3654,]    7
[3655,]    8
[3656,]    9
[3657,]    8
[3658,]    8
[3659,]    6
[3660,]    7
[3661,]    8
[3662,]    5
[3663,]    7
[3664,]    3
[3665,]    4
[3666,]    7
[3667,]    9
[3668,]    6
[3669,]    7
[3670,]    5
[3671,]    5
[3672,]    7
[3673,]    6
[3674,]    7
[3675,]    5
[3676,]    6
[3677,]    4
[3678,]    5
[3679,]    3
[3680,]    4
[3681,]    8
[3682,]    6
[3683,]    7
[3684,]    5
[3685,]    3
[3686,]    7
[3687,]    9
[3688,]    5
[3689,]    5
[3690,]    5
[3691,]    6
[3692,]    6
[3693,]    7
[3694,]    7
[3695,]    9
[3696,]    6
[3697,]    4
[3698,]    4
[3699,]    4
[3700,]    6
[3701,]    4
[3702,]    3
[3703,]    5
[3704,]    3
[3705,]    3
[3706,]    3
[3707,]    9
[3708,]    3
[3709,]    8
[3710,]    4
[3711,]    6
[3712,]    4
[3713,]    3
[3714,]    5
[3715,]    8
[3716,]    7
[3717,]    9
[3718,]    8
[3719,]    7
[3720,]    6
[3721,]    7
[3722,]    5
[3723,]    9
[3724,]    8
[3725,]    5
[3726,]    7
[3727,]    6
[3728,]    3
[3729,]    5
[3730,]    5
[3731,]    7
[3732,]    6
[3733,]    5
[3734,]    8
[3735,]    7
[3736,]    4
[3737,]    4
[3738,]    8
[3739,]    7
[3740,]    4
[3741,]    5
[3742,]    5
[3743,]    5
[3744,]    4
[3745,]    3
[3746,]    4
[3747,]    6
[3748,]    8
[3749,]    4
[3750,]    8
[3751,]    7
[3752,]    5
[3753,]    6
[3754,]    4
[3755,]    6
[3756,]    7
[3757,]    8
[3758,]    8
[3759,]    7
[3760,]    5
[3761,]    6
[3762,]    7
[3763,]    7
[3764,]    7
[3765,]    5
[3766,]    3
[3767,]    6
[3768,]    7
[3769,]    1
[3770,]    5
[3771,]    6
[3772,]    6
[3773,]    4
[3774,]    5
[3775,]    7
[3776,]    7
[3777,]    7
[3778,]    6
[3779,]    6
[3780,]    6
[3781,]    6
[3782,]    3
[3783,]    6
[3784,]    7
[3785,]    3
[3786,]    3
[3787,]    3
[3788,]    4
[3789,]    4
[3790,]    7
[3791,]    4
[3792,]    2
[3793,]    8
[3794,]    7
[3795,]    5
[3796,]    6
[3797,]    8
[3798,]    6
[3799,]    7
[3800,]    4
[3801,]    9
[3802,]    7
[3803,]    7
[3804,]    9
[3805,]    5
[3806,]    2
[3807,]    4
[3808,]    3
[3809,]    4
[3810,]    1
[3811,]    5
[3812,]    7
[3813,]    3
[3814,]    7
[3815,]    6
[3816,]    3
[3817,]    3
[3818,]    6
[3819,]    6
[3820,]    3
[3821,]    8
[3822,]    7
[3823,]    2
[3824,]    7
[3825,]    8
[3826,]    4
[3827,]    3
[3828,]    6
[3829,]    6
[3830,]    5
[3831,]    4
[3832,]    6
[3833,]    3
[3834,]    5
[3835,]    4
[3836,]    5
[3837,]    7
[3838,]    6
[3839,]    6
[3840,]    5
[3841,]    3
[3842,]    6
[3843,]    5
[3844,]    4
[3845,]    5
[3846,]    8
[3847,]    7
[3848,]    8
[3849,]    8
[3850,]    5
[3851,]    4
[3852,]    5
[3853,]    6
[3854,]    5
[3855,]    9
[3856,]    6
[3857,]    7
[3858,]    7
[3859,]    9
[3860,]    8
[3861,]    8
[3862,]    8
[3863,]    8
[3864,]    6
[3865,]    8
[3866,]    7
[3867,]    7
[3868,]    4
[3869,]    6
[3870,]    6
[3871,]    8
[3872,]    7
[3873,]    3
[3874,]    8
[3875,]    6
[3876,]    5
[3877,]    5
[3878,]    7
[3879,]    6
[3880,]    8
[3881,]    6
[3882,]    5
[3883,]    3
[3884,]    5
[3885,]    5
[3886,]    8
[3887,]    6
[3888,]    8
[3889,]    7
[3890,]    6
[3891,]    6
[3892,]    7
[3893,]    8
[3894,]    7
[3895,]    7
[3896,]    5
[3897,]    4
[3898,]    3
[3899,]    5
[3900,]    9
[3901,]    9
[3902,]    8
[3903,]    7
[3904,]    7
[3905,]    7
[3906,]    5
[3907,]    4
[3908,]    6
[3909,]    6
[3910,]    5
[3911,]    7
[3912,]    6
[3913,]    4
[3914,]    6
[3915,]    6
[3916,]    8
[3917,]    6
[3918,]    6
[3919,]    8
[3920,]    4
[3921,]    4
[3922,]    5
[3923,]    8
[3924,]    8
[3925,]    7
[3926,]    7
[3927,]    4
[3928,]    6
[3929,]    6
[3930,]    6
[3931,]    6
[3932,]    7
[3933,]    4
[3934,]    6
[3935,]    6
[3936,]    7
[3937,]    3
[3938,]    8
[3939,]    4
[3940,]    4
[3941,]    7
[3942,]    6
[3943,]    4
[3944,]    6
[3945,]    6
[3946,]    6
[3947,]    6
[3948,]    6
[3949,]    4
[3950,]    7
[3951,]    7
[3952,]    7
[3953,]    3
[3954,]    5
[3955,]    2
[3956,]    3
[3957,]    4
[3958,]    6
[3959,]    2
[3960,]    7
[3961,]    5
[3962,]    6
[3963,]    6
[3964,]    6
[3965,]    3
[3966,]    4
[3967,]    7
[3968,]    8
[3969,]    5
[3970,]    6
[3971,]    5
[3972,]    5
[3973,]    5
[3974,]    6
[3975,]    4
[3976,]    8
[3977,]    5
[3978,]    2
[3979,]    3
[3980,]    7
[3981,]    6
[3982,]    5
[3983,]    7
[3984,]    6
[3985,]    4
[3986,]    6
[3987,]    7
[3988,]    9
[3989,]    5
[3990,]    7
[3991,]    4
[3992,]    3
[3993,]    2
[3994,]    5
[3995,]    6
[3996,]    6
[3997,]    8
[3998,]    9
[3999,]    6
[4000,]    6
data.frame(obs = ppd) %>% 
  ggplot(aes(x=obs)) +
  geom_histogram() +
  labs(x="Observed water (out of 9)")

Sampling parameters from the prior

Oops, we’ve jumped ahead of ourselves! Best practice is to simulate values from your prior first and check to see if those priors are reasonable.

m1p <-
  brm(data = list(w = 6),                            
      family = binomial(link = "identity"),          
      w | trials(9) ~ 0 + Intercept,                 
      prior(beta(1, 1), class = b, lb = 0, ub = 1),  
      iter = 5000, warmup = 1000, seed = 3, chains=1,          
      sample_prior = "only")
samples_from_prior = as_draws_df(m1p)
samples_from_prior %>% 
  ggplot(aes(x=b_Intercept)) +
  geom_density(fill = "grey", color = "white") +
  labs(x="Proportion water", title="Prior")

Simulating observations from the prior

We may also want the PRIOR PREDICTIVE DISTRIBUTION which is the expected observiations given our prior.

prior_pd = posterior_predict(m1p)
data.frame(obs = prior_pd) %>% 
  ggplot(aes(x=obs)) +
  geom_histogram() +
  labs(x="Observed water (out of 9)")

Simulating from your priors – prior predictive simulation – is an essential part of modeling. This allows you to see what your choices imply about the data. You’ll be able to diagnose bad choices.

an aside about learning in R

At this point in the course, I’m going to start throwing a lot of code at you. Do I expect you to memorize this code? Of course not.

Do you need to understand every single thing that’s happening in the code? Nope.

But, you’ll learn a lot by taking the time to figure out what’s happening in a code chunk. Class time will frequently include exercises where I ask you to adapt code I’ve shared in the slides to a new dataset or to answer a new problem. When doing so, go back through the old code and figure out what’s going on. Run the code one line at a time. Always observe the output and take some time to look at the object that was created or modified. Here are some functions that will be extremely useful:

str() # what kind of object is this? what is its structure?
dim() # what are the dimensions (rows/columns) of this object
head() # give me the first bit of this object
str(prior_pd)
 int [1:4000, 1] 6 1 1 0 9 8 0 2 4 3 ...
 - attr(*, "dimnames")=List of 2
  ..$ : NULL
  ..$ : NULL
dim(prior_pd)
[1] 4000    1
head(prior_pd)
     [,1]
[1,]    6
[2,]    1
[3,]    1
[4,]    0
[5,]    9
[6,]    8

Continous outcomes

The globe tossing example is cute and easy to work with, but let’s move towards the kinds of variables we more frequently work with. Let’s create a model for some outcome, \(y\) that is continuous.

\[\begin{align*} y_i &\sim \text{Normal}(\mu, \sigma) \\ \mu &\sim \text{Normal}(0, 10) \\ \sigma &\sim \text{Uniform}(0, 5) \end{align*}\]

set.seed(9)
y = rnorm(n = 31, mean = 4, sd = .5)
m2 = brm(
  data = list(y=y),
  family = gaussian,
  y ~ 1,
  prior = c(prior( normal(0,10), class=Intercept),
            prior( uniform(0,5), class=sigma)),  
      iter = 5000, warmup = 1000, seed = 3, chains=1,
  file = here("files/models/m21.2")
)

An example: weight and height

Using the Howell data (don’t load the rethinking package because it can interfere with brms).

data("Howell1", package = "rethinking")
d <- Howell1
str(d)
'data.frame':   544 obs. of  4 variables:
 $ height: num  152 140 137 157 145 ...
 $ weight: num  47.8 36.5 31.9 53 41.3 ...
 $ age   : num  63 63 65 41 51 35 32 27 19 54 ...
 $ male  : int  1 0 0 1 0 1 0 1 0 1 ...
library(measurements)
d$height <- conv_unit(d$height, from = "cm", to = "feet")
d$weight <- conv_unit(d$weight, from = "kg", to = "lbs")
rethinking::precis(d)
             mean         sd      5.5%    94.5%      histogram
height  4.5362072  0.9055921  2.661042   5.4375      ▁▁▁▁▂▂▇▇▁
weight 78.5079631 32.4502290 20.636856 120.1583 ▁▂▃▂▂▁▁▃▅▇▇▃▂▁
age    29.3443934 20.7468882  1.000000  66.1350      ▇▅▅▃▅▂▂▁▁
male    0.4724265  0.4996986  0.000000   1.0000     ▇▁▁▁▁▁▁▁▁▇
d2 <- d[ d$age >= 18, ]

exercise

Write a mathematical model for the weights in this data set. (Don’t worry about predicting from other variables yet.)

solution

\[\begin{align*} w &\sim \text{Normal}(\mu, \sigma) \\ \mu &\sim \text{Normal}(130, 20) \\ \sigma &\sim \text{Uniform}(0, 25) \\ \end{align*}\]

\[\begin{align*} w &\sim \text{Normal}(\mu, \sigma) \\ \mu &\sim \text{Normal}(130, 20) \\ \sigma &\sim \text{Uniform}(0, 25) \\ \end{align*}\]

exercise

Simulate from your priors (parameters values and prior predictive values).

solution

Sample from your priors:

m3p = brm(
  data = d2,
  family = gaussian,
  weight ~ 1,
  prior = c(prior( normal(130,20), class=Intercept),
            prior( uniform(0,25), class=sigma, lb=0, ub=25)),  
      iter = 5000, warmup = 1000, seed = 3, chains=1,
  sample_prior = "only")

Sampling parameter estimates for the Intercept.

pairs(m3p)

This is a different (shorter) way to plot your posterior. Good things: it automatically includes the scatterplot so you can see the implications of how these parameters correlate. Bad things: not customizable and not useable when you have a lot of parameters.

Simulate values of weight.

prior_pd = posterior_predict(m3p)
dim(prior_pd)
[1] 4000  352
as.data.frame(prior_pd) %>% 
  pivot_longer(everything()) %>% 
  ggplot(aes(x=value)) +
  geom_histogram() +
  labs(x="Expected observed weights (based on prior)")

Another shorter way:

pp_check(m3p)

Fit the model

m3 = brm(
  data = d2,
  family = gaussian,
  weight ~ 1,
  prior = c(prior( normal(130,20), class=Intercept),
            prior( uniform(0,25), class=sigma, lb=0, ub=25)),  
      iter = 5000, warmup = 1000, seed = 3, chains=1,
  file = here("files/models/m21.3"))
posterior_summary(m3)
                Estimate  Est.Error         Q2.5        Q97.5
b_Intercept    99.221909 0.75720356    97.751074   100.737488
sigma          14.291987 0.55198322    13.254363    15.428183
Intercept      99.221909 0.75720356    97.751074   100.737488
lprior         -8.318377 0.05827127    -8.433538    -8.203915
lp__        -1441.294276 1.02695551 -1444.235605 -1440.292326
pairs(m3)
Code
posterior_predict(m3) %>% 
  as.data.frame() %>% 
  pivot_longer(everything()) %>% 
  ggplot(aes(x=value)) +
  geom_density(fill = "grey", color = "white") +
  geom_density( aes(x = weight), data=d2, inherit.aes = F) 
pp_check(m3)

Adding in a linear component

We might assume that height and weight are associated with each other. Indeed, within our sample:

plot(d2$weight ~ d2$height)

exercise

Update your mathematical model to incorporate height.

\[\begin{align*} w_i &\sim \text{Normal}(\mu_i, \sigma) \\ \mu_i &= \alpha + \beta h_i \\ \alpha &\sim \text{Normal}(??, ??) \\ \beta &\sim \text{Normal}(0, 25) \\ \sigma &\sim \text{Uniform}(0, 25) \\ \end{align*}\]

exercise

Update your mathematical model to incorporate height.

\[\begin{align*} w_i &\sim \text{Normal}(\mu_i, \sigma) \\ \mu_i &= \alpha + \beta (h_i - \bar{h}) \\ \alpha &\sim \text{Normal}(130, 20) \\ \beta &\sim \text{Normal}(0, 25) \\ \sigma &\sim \text{Uniform}(0, 25) \\ \end{align*}\]

\[\begin{align*} w_i &\sim \text{Normal}(\mu_i, \sigma) \\ \mu_i &= \alpha + \beta (h_i - \bar{h}) \\ \alpha &\sim \text{Normal}(130, 20) \\ \beta &\sim \text{Normal}(0, 25) \\ \sigma &\sim \text{Uniform}(0, 25) \\ \end{align*}\]

To update our brms code:

d2$height_c = d2$height - mean(d2$height)
  
m4p = brm(
  data = d2,
  family = gaussian,
  weight ~ 1 + height_c,
  prior = c(prior( normal(130,20), class=Intercept),
            prior( normal(0,25),   class=b),
            prior( uniform(0,25),  class=sigma, lb=0, ub=25)),  
      iter = 5000, warmup = 1000, seed = 3, chains=1,
  sample_prior = "only")
set.seed(9)
samples_from_prior = as_draws_df(m4p)
str(samples_from_prior)
draws_df [4,000 × 9] (S3: draws_df/draws/tbl_df/tbl/data.frame)
 $ b_Intercept: num [1:4000] 167 155.6 149.6 88.7 176.4 ...
 $ b_height_c : num [1:4000] -57.9 -57.2 -18.5 52 -47.3 ...
 $ sigma      : num [1:4000] 14.1 11.8 13.3 2.1 21.1 ...
 $ Intercept  : num [1:4000] 167 155.6 149.6 88.7 176.4 ...
 $ lprior     : num [1:4000] -15.7 -14.7 -12 -15.6 -15.8 ...
 $ lp__       : num [1:4000] -13.8 -12.9 -10.2 -14.9 -14.6 ...
 $ .chain     : int [1:4000] 1 1 1 1 1 1 1 1 1 1 ...
 $ .iteration : int [1:4000] 1 2 3 4 5 6 7 8 9 10 ...
 $ .draw      : int [1:4000] 1 2 3 4 5 6 7 8 9 10 ...
Code
d2 %>% 
  ggplot(aes(x=height_c, y=weight)) +
  geom_blank() +
  geom_abline( aes(intercept=b_Intercept, slope=b_height_c), 
               data=samples_from_prior[1:50, ],
               alpha=.3) +
  scale_x_continuous(name = "height(feet)", 
                     breaks=seq(4,6,by=.5)-mean(d2$height),
                     labels=seq(4,6,by=.5))

Describe in words what’s wrong with our priors.

Slope should not be negative. How can we fix this?

Could use a uniform distribution bounded by 0.

m4p = brm(
  data = d2,
  family = gaussian,
  weight ~ 1 + height_c,
  prior = c(prior( normal(130,20), class=Intercept),
            prior( uniform(0,25),   class=b),
            prior( uniform(0,25),  class=sigma, lb=0, ub=25)),  
      iter = 5000, warmup = 1000, seed = 3, chains=1,
  sample_prior = "only")

exercise

Fit the new weight model to the data.

solution

m4 = brm(
  data = d2,
  family = gaussian,
  weight ~ 1 + height_c,
  prior = c(prior( normal(130,20), class=Intercept),
            prior( uniform(0,25),   class=b),
            prior( uniform(0,25),  class=sigma, lb=0, ub=25)),  
      iter = 5000, warmup = 1000, seed = 3, chains=1,
  file=here("files/models/m21.4"))
posterior_summary(m4)
               Estimate Est.Error        Q2.5       Q97.5
b_Intercept    99.21508 0.5154139    98.17317   100.20955
b_height_c     24.74661 0.2609169    24.05150    24.99355
sigma          10.36090 0.3898531     9.65234    11.12228
Intercept      99.21508 0.5154139    98.17317   100.20955
lprior        -11.53739 0.0396881   -11.61861   -11.46176
lp__        -1332.15882 1.3936194 -1335.90639 -1330.52002

exercise

Draw lines from the posterior distribution and plot with the data.

solution

Code
set.seed(9)
samples_from_posterior = as_draws_df(m4)
d2%>% 
  ggplot(aes(x=height_c, y=weight)) +
  geom_point(size=.5) +
  geom_abline( aes(intercept=b_Intercept, slope=b_height_c), 
               data=samples_from_posterior[1:50, ],
               alpha=.3,
               color="#1c5253") +
  scale_x_continuous(name = "height(feet)", 
                     breaks=seq(4,6,by=.5)-mean(d2$height),
                     labels=seq(4,6,by=.5))

A side note: a major concern or critique of Bayesian analysis is that the subjectivity of the priors allow for nefarious behavior. “Putting our thumbs on the scale,” so to speak. But priors are quickly overwhelmed by data. Case in point:

m4e = brm(
  data = d2,
  family = gaussian,
  weight ~ 1 + height_c,
  prior = c(prior( normal(130,20), class=Intercept),
            prior( normal(-5,5),   class=b),
            prior( uniform(0,25),  class=sigma, lb=0, ub=25)),  
      iter = 5000, warmup = 1000, seed = 3, chains=1,
  file=here("files/models/m21.4e"))
posterior_summary(m4e)
                Estimate Est.Error         Q2.5       Q97.5
b_Intercept    99.197733 0.5035653    98.214092   100.15932
b_height_c     35.775818 1.8832462    32.177572    39.50070
sigma           9.529429 0.3687178     8.839379    10.27233
Intercept      99.197733 0.5035653    98.214092   100.15932
lprior        -44.172476 3.0738968   -50.456409   -38.49994
lp__        -1334.629214 1.1849503 -1337.645046 -1333.23509

You’ll only really get into trouble with uniform priors that have a boundary, if true population parameter is outside your boundary. A good rule of thumb is to avoid the uniform distribution. We’ll cover other options for priors for \(\sigma\) in future lectures, but as a preview, the exponential distribution works very well for this!